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lO Abstract 

A perturbation framework, called the quasi-stationary analysis (QSA), 
I I is developed to analyze metastable behavior in stochastic processes with 

Qi^ random internal and external states. The QSA is illustrated with a model 

•^r of gene expression that displays bistable switching. In this model, the 

external state represents the number of protein molecules produced by 

a hypothetical gene. Once produced, a protein is eventually degraded. 

C^ The internal state represents the activated or unactivated state of the 

C gene; in the activated state the gene produces protein more rapidly than 

I the unactivated state. The gene is activated by a dimer of the protein 

it produces so that the activation rate depends on the current protein 
^ level. This is a well studied model, and several model reductions and 

^ diffusion approximation methods are available to analyze its behavior. 

However, it is unclear if these methods accurately approximate long-time 
metastable behavior (i.e., mean switching time between metastable states 
of the bistable system). Diffusion approximations are generally known to 
fail in this regard. It is shown that a diffusion approximation based on 
\l a quasi-steady-state reduction (stochastic averaging), which averages the 

internal state out and reduces the process to a continuous Markov process 



^ 



in 



o 



_ ■ for the external state, provides unreliable accuracy. On the other hand, 

. . the QSA approximation is consistently accurate. 



X 



1 Introduction 

A common feature found in many stochastic models of biological processes is 
a distinction between internal and external states [36]. There are numerous 
examples of such Markov processes used as models for biological phenomena 
[2, 10, 17, 27, 31]. Examples of an internal state include the number of open ion 
channels in the membrane of a neuron that affect its membrane voltage [16] 
and the on/off state of a gene that affects its protein production rate ]28] . The 
distinction between internal and external states should not be confused with the 
concept of 'intrinsic' and 'extrinsic' noise (see Ref. [35[ for an example related 
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to gene expression). A system with internal degrees of freedom is a classical 
idea in physics and applied mathematics, and the extension of this concept to 
Markov processes with internal states is well known in the literature [14, 18,20]. 
The canonical example is the velocity jump process, which has a continuous 
external state, usually representing the physical location of a particle or some 
other observable quantity, that evolves deterministically with a velocity that 
depends on the current value of the internal state. The internal state evolves 
randomly according to a discrete jump Markov process. The simplest velocity 
jump process is called the telegraph process. The external state X{t) evolves 
according to the stochastic differential equation 

X{t)=X{0)+ [ v{S{T))dT, (1.1) 

Jo 

where the velocity w(— 1) — — w- and v{l) = v+. The internal state S{t) can 
take the value ±1 at any given time and is governed by the discrete Markov 
process 

S{t) = S{0) + Y,{[ ax{S{T))dT) - Y2{ f /3(1 - x(5(T)))dr), (1.2) 

Jo Jo 

where Yi2{t) are independent, unit Poisson processes, and the indicator function 
x(s) = 1 if s = 1 and is zero otherwise. The transition rates (or intensities) 
a and (5 determine how quickly the particle switches velocity. Velocity jump 
processes have been used to model molecular motor transport [29], dynamic 
microtubule polymerization ]2], and chemotaxis ]31]. 

The Markov process becomes nonlinear when internal state transitions de- 
pend on the external state. For example, the process {X{t), S{t)) is nonlinear if 
the transition rates depend on X{t) (i.e., a — ^ a{X{t)) and /? —> f3{X{t))). In- 
terestingly, nonlinearities can lead to metastable behavior under so-called weak 
noise conditions. For the above example, this occurs when the transition rates 
are large compared to v/L, where L is a characteristic space scale on which 
the process is observed. In the large transition rate limit, the process becomes 
deterministic and 

a{x) \ / (3{x) \ 

a(.T)+/3(a;)j^'+" [a{x) + f3{x) ) ""-- ^^^ ^ 

Metastable behavior is important because it represents fluctuation-induced phe- 
nomena that cannot be predicted by the deterministic limit. The standard ex- 
ample of metastable behavior is Brownian motion in a double well potential. On 
short timescales the particle is most likely found near one of the two minima, 
and on long timescales the particle can transition over the energy barrier that 
separates the two wells. Other types of Markov processes, such as the velocity- 
jump process, can behave like diffusion in a double well potential. Here, we 
consider processes with an internal state that evolves according to a discrete 
jump Markov process. But unlike the example velocity jump process described 
above, we also consider noise in the external state. 



Assume that there is a well-defined limit where both noise sources are re- 
moved and the process converges to a deterministic system describing the ex- 
ternal state. Under weak noise conditions, the dynamics of the deterministic 
system strongly influence the dynamics of the stochastic process. In particular, 
we are interested in the case where the deterministic system has multiple stable 
solutions depending on the initial conditions. On short timescales, a trajec- 
tory of the stochastic process fluctuates about the deterministic trajectory that 
has the same initial conditions. However, metastable behavior in the stochastic 
process is not seen in the deterministic system because it depends on a small 
amount of noise present in the system to cause a transition from one of the sta- 
ble deterministic solutions to another. Metastable transitions occur on a long 
timescale. 

Metastable transitions by nonlinear Markov processes with an internal and 
external state are more difficult to analyze than diffusion in a potential well, and 
exact analytical solutions are rarely possible. Monte Carlo simulation does not 
work either because metastable behavior requires far to much processor time to 
be practical. It is therefore necessary to develop approximation methods. One 
way to approximate is to eliminate a noise source. To eliminate noise in the 
internal state, an adiabatic limit is taken where the transitions are infinitely 
rapid so that the internal state can be averaged out. Eliminating noise in the 
external state results in a velocity jump process where the external state evolves 
deterministically in between random jumps in the internal state. However, the 
timescale for metastable transitions is very sensitive to both the type of noise 
and the noise strength, and eliminating a noise source can lead to large errors. 

Another way to reduce the complexity of the model is with a quasi-steady- 
state (QSS) reduction. This is very similar to the adiabatic limit, but uses a 
perturbation approach so that higher order terms can be included that account 
for noise in the internal state. The QSS reduction includes effects from both 
noise sources, but approximates the coupled (internal, external) process with 
a single continuous Markov process for the external state. This approximation 
reduces the problem to diffusion in a double potential well. The underlying 
assumption behind the approximation is that the random internal state is well 
approximated by a random variable chosen from its steady-state distribution 
conditioned on a fixed external state. Intuitively, the assumption is valid if the 
internal state fluctuates much more rapidly than the external state (i.e., the 
small noise conditions described above) . 

While the QSS reduction is a useful tool in many circumstances, it is not 
accurate for characterizing metastability. Recently, in a model of molecular mo- 
tor transport, metastable behavior was observed in a velocity jump process [27]. 
Following this observation, alternative perturbation methods were developed to 
analyze the metastable behavior [16,28,30]. However, these methods only work 
when the external state evolves deterministically as in the example described 
above. In this paper, we further develop this theory to account for noise in the 
external state. 

To describe metastable behavior, it is necessary to approximate both the 
effective potential for the process and various timescales for metastable transi- 



tions. For ID continuous Markov processes, quantities such as the potential are 
straightforward to define, and if a stationary solution exists it must have zero 
probability flux. Given its usefulness at describing the qualitative features, we 
would like to know if we can define an effective potential in general. For higher 
dimensional continuous Markov processes, the potential well is no longer well 
defined when the drift velocity field is not the gradient of a potential. This is 
closely related to detailed balance conditions and thermodynamic equilibrium. 
In higher dimensions, it is possible for the stationary density to exhibit a nonzero 
probability flux. Developing a systematic formalism to describe nonequilibrium 
stationary behavior is particularly relevant in biology. It turns out that an effec- 
tive potential can still be defined through perturbation theory [32] and transition 
path theory [13,25]. These tools can also be used to approximate the timescale 
associated with metastable transitions. The methods presented here fit within 
the perturbation framework. 

The theory of large deviations [8,9,33[ is the mathematical foundation for the 
techniques used to study metastable transitions (rare events). Here, we focus 
on perturbation-theory-based techniques [32[, which we refer to as the quasi- 
stationary analysis (QSA). The QSA was developed to analyze the differential 
Kolmogorov equation [11[, which describes the process through its probability 
density function. For a continuous Markov process, the QSA is well-developed 
for the Fokker-Planck equation [22-24, 26, 32, 34[. The QSA has also be apphed 
to the Master equation to analyze certain birth-death processes [1,3-7, 12, 15, 37[. 
(The Fokker-Planck and Master equation are instances of the more general 
differential Kolmogorov equation.) However, as far as we are aware, no one has 
applied these methods to study weak noise problems where QSS reduction (i.e., 
stochastic averaging) and the large system size limit are both necessary to reach 
the deterministic system. Furthermore, large deviation theory does not provide 
a means of explicitly calculating the preexponential factor (see Section 4.1), 
which is a significant part of the leading order transition time and stationary 
density approximations. 

There are several advantages to the QSA over diffusion-approximation and 
model-reduction techniques. First, if the process includes a discrete state, the 
QSA provides an approximation that accounts for all moments of the jump prop- 
agator, whereas the diffusion approximations include only the first two moments 
(e.g., a diffusion approximation of a discrete jump Markov process by trunca- 
tion of a Kramers-Moyal (KM) expansion). Second, the approximation provides 
a uniformly accurate approximation of the stationary probability density func- 
tion. Third, physically meaningful quantities, such as the effective potential and 
metastable transition rates, can be generalized to processes that do not assume 
detailed balance. Finally, the QSA can be applied to high dimensional (in the 
external state) problems displaying nonequilibrium stationary behavior. 

In this paper, we develop the QSA for a general class of Markov processes 
that have a discrete internal state, and we illustrate the analysis using a simple 
model problem. The model problem is ideal because it is possible to derive a 
QSS diffusion approximation and several different model reductions, allowing us 
to (i) develop the general QSA procedure for different types of Markov processes 



(ii) compare the resulting approximations with the QSS diffusion approximation. 
In particular, we are interested in approximating two quantities: the timescales 
for metastable transitions and the effective potential. The analysis of the model 
problem should inform our understanding about when reduction techniques, 
such as a diffusion approximation, fail to approximate these two quantities and 
why. Previous work has shown that diffusion approximations lead to errors 
in both the adiabatic [12,38] and quasi-deterministic limits [28]. But what 
happens when both noise sources are present? When is one noise source more 
significant than the other? If the QSS reduction fails, why does it fail? Is it 
due to large deviation errors like the system-size expansion, or is it because the 
QSS assumption is invalid? Does it fail for the same reasons in each limit? 

The paper is organized as follows. First, in Section 2 we introduce the model 
problem along with various approximations and model reductions. In Section 3 
we describe the dynamics of the deterministic limit and explore the two limits 
where noise in the either the internal or external state is removed. The QSA is 
presented in Section 4, and results are shown in Section 5. 

2 Model problem 

Consider the following as an example of a discrete Markov process with an 
internal state. A population of proteins changes stochastically according to a 
jump Markov process, where the protein production rate depends on the internal 
state of the system. The hypothetical gene responsible for producing the protein 
has a promotor that can be turned on or off. When the promotor is on, protein 
is produced at higher rate than when it is off. For simplicity all parameters are 
presented in nondimensional form. The following state diagram, where iV„ is 
the state where n proteins are present in the system, represents the external 
state transitions: 

T(S(t)) T{S{t)) x(S(t)) T(S(t)) 

N^ ^ Ni ^ N2--- ^ Nr, ^ ■■■ , (2.1) 

6 26 n5 (n + l)i5 

where we set 5 — 1 and S{t) is the two state stochastic process governing the 
on/off state of the promotor; S{t) — 1 when the promoter is on and S{t) — 
when it is off. The production rate is a function of the promotor state, with 

t(0) = aae, t(1) = a,. (2.2) 

The nondimensional parameter a controls how much spontaneous protein pro- 
duction occurs when the promotor is off, and we assume < cr < 1 so that 
protein production is higher when the promotor is on. To get nonlinear phe- 
nomena, the internal state transitions must depend on the external state. A 
simple model of the promotor fluctuations is given by 

(off) ;i^ (on), (2.3) 

Qi/3 



where N{t) is the number of protein copies. This is a model of a promotor that 
can be turned on by a dimer of the protein product [17]. 

To fully model the dimerization reaction, an additional state variable is re- 
quired which expands the state space of the process. Under the assumption that 
there are much fewer dimers than monomers, a quasi-steady-state approxima- 
tion can be made [17, 19]. Of course, this is a model reduction approximation, 
and we do not fail to recognize that this is a paper about the validity model 
reduction techniques. However, because we are interested in a simple model 
we can use to explore various mathematical issues, we accept (2.3) as the basic 
model of promotor switching, though we make no claims about its accuracy in 
approximating the dimerization reaction. 

The transitions between the two promotor states are assumed to be fast by 
specifying that l/ai <C 1 is a small parameter. Thus, the full Markov process 
for the combined external/internal state of the system is M{t) = {N{t),S{t)). 
Note that neither process is Markovian by itself. The Chapman-Kolmogorov 
(CK) equation governing the probability density function pj(n,t) = p{j^n,t \ 
jo, no, to) for M{t) is 

d f n^ \ 

— Po(n,t) = [(E+ - l)n + aofT(E- - 1)] po + a; f -— po -l- /3pi j (2.4a) 

-pi(n, t) = [(E+ - l)n + ae(E- - 1)] Pi -|- ai f ^Po - /3pi j , (2.4b) 

where the jump operators E are defined by 

E±/(r^)-/(n±l). (2.5) 

It is convenient to introduce the following notation. We can rewrite the 
probability density function as 

Prob[5(i) = s, X{t) C,{x,x + dx)] = Prob[5(t) = s\X{t) e (x, a; + dx)] 

X Prob[X(t) e {x,x + dx)]. (2.6) 
Define the conditional internal state distribution to be 

w,{x,t) = Vvoh[s = S{t)\X{t) e {x,x + dx)i (2.7) 

and the marginal external state density to be 

u{x, t) = Vroh[X{t) e{x,x + dx)]. (2.8) 

Introduce vector notation for the probability density with 

p(a:,t) = (po(a:,i),pi(x,t))^. (2.9) 

We write the stationary versions of (2.7)-(2.9) as 

lim p(a;,i) = p(x) = w(a;)u(x), (2-10) 

where w(a;) is a 2-vector and 

Y,Ws{x) = l. (2.11) 



2.1 Semi-continuous process 

A closely-related process assumes an external state that evolves according to a 
continuous Markov process instead of a discrete birth-death process as in the 
previous section. Consider the Langevin equation 

dX{t) = -v{S{t))dt + y/b{S{t))dW{t). (2.12) 

where S{t) is again the internal state. While this process can be quite general, 
in order to make the connection to the discrete process (2.1) we define it as 
the diffusion approximation via a truncated KM expansion. The mean number 
of proteins when the promotor is on is ac. When a^ is large we can rescale 
to a continuous variable X{t) = N{t)/ae. The internal state S{t) then evolves 
according to 

(off) "^ (on), (2.13) 



It is straight forward to show that the drift in each state is 

i;(0)=(T-X, v{l)^l-X, (2.14) 

and the diffusivities are 

6(0) = — ((7 + X), 6(1) = — (1 + X). (2.15) 

The full Markov process is then C{t) = {X{t)^ S{t)), and its corresponding dif- 
ferential Kolmogorov equation is 

d d I d"^ 

—po{x,t) = - — [(ct-x)po] + ^^[(^ + ^)Po] -"i (^^^Po - Ppi) 

(2.16a) 

d d 1 (9^ 

— Pi(a;,i) = - — [(l-a;)Pi] + ^ — q^[(1 + ^)pi] + "i (a;^Po - ^Pi) ■ 



(2.16b) 



2.2 Diffusion approximation 



There is a standard asymptotic method available to analyze both processes 
introduced above, based on the same principle as that of a quasi-steady-state 
(QSS) reduction. If the system transitions between the internal states much 
more quickly than it transitions between its external states, then the internal 
states are approximately given by their steady state distribution conditioned on 
a fixed external state. The combined process {X{t), S{t)), can be reduced using 
a projection method to eliminate (average out) the internal states S(t) to obtain 
a Markov process that approximates the external state X(t). In other words, 
although X{t) is not strictly Markovian, it may be approximately Markovian. 



A projection method results in a scalar Fokker-Planck equation for the 
marginal external-state probability density function 



l{x,t) =^ps{x,t). 



(2.17) 



s=0 



For a general discussion of the QSS projection method see [11]. For brevity we 
only quote the result here; for further details of the QSS reduction of the model 
problem considered here, we refer the reader to [17]. The result is 



du 



l((vW> ") + ^^ ((b(.)> .) + ^1: (Di.)'£) , (2.18) 



where 



(v(a;)) 



I3{a-x) x^{l-x) 



p + x^ 



(3 + x^ 



a; f /3{a + x) , x^{l + x) 



D{x) 



P{{a - x) - {^(x))){a - x) , x\{\-x)-{^{x))){\-x) 



(/3 + 2;2) 



2^2 



(/? + x2) 



2\2 



(2.19) 
(2.20) 

(2.21) 



Protein fluctuations are captured by (}o{x)) and proniotor fluctuations by D{x). 
We can rewrite the Fokker-Planck equation (2.18) in conservation form, with 



du 



dj 

dx ' 



J{x,t) = a{x)u{x,t) — B{x) 



du 
dx' 



where 



a{x) = (v(:r)) - ^'^^, B{x) = - ((b(x)) + D{x)) . 



(2.22) 



(2.23) 



The mean first passage time (MFPT) for any given starting position xo to reach 
x± is given exactly by 



T = 



^± ry pe(z)-efe) 



xq Jo 



B{v) 



-dzdy, Q{x) 



ajy) 

B{y) 



dy. 



(2.24) 



To make comparisons to other approximations, we make the following defini- 
tions. Let 



$(a;) 



Hy)) 



,^ My)) + D{y) '^^' ^^"^^ J,^ (b(y)) + D{y) 



rJHy)) 



dy 



dy. 



The MFPT is then approximated by T± = T{x±) ^ 1/A±, where 
^^ ^ (B{x^^^^„^^^^^^,^^ e-*(-*)exp \--^x.) 



(2.25) 



(2.26) 



3 Limiting processes 

In the previous section, we define the model problem. We call the full version 
presented in Section 2 the discrete process, which makes no assumption about 
the size of the parameters ai or ac. The semi-continuous process from Section 
2.1 is an approximation of the discrete process where x — n/ac is treated as 
a continuous variable and we assume that ac ^ 1 is a large parameter; no 
assumption is made about the size of a;. If we assume that ai is also a large 
parameter then further reduction is possible using a diffusion approximation, 
call it the QSS process, which is presented in Section 2.2. All three versions 
contain terms that depend on a; and a^, and if these parameters are assumed 
to be large, all three should account for contributions of noise in the external 
and internal state. 

Of course, further model reduction is possible by removing one source of 
noise: either ai — > oo or ac — >■ oo. The former is known in the literature as the 
adiabatic limit (see [17,38]), and the later we call the quasi-deterministic (QD) 
limit. If both limits are taken, a well defined deterministic dynamical system 
is obtained. Note that all three versions of the model problem — the discrete, 
semi-continuous, and QSS processes — converge to the same deterministic limit 
(3.1). In the rest of this section we explore each limit in turn. 

3.1 Deterministic limit ai — )■ oo and «c — > oo 

If we take the limit ac -^ oo and ofj — > oo, the resulting deterministic dynamics 
is 

. _ P{a-x) + x\l-x) 

P + x^ ■ ^•^■^> 

Assuming that cr ^ 1, the system is described as follows. For /3_ < /3 < /3_|_, 
where 

P^^4a + 0ia'), /3+^i + |+0(a2) (3.2) 

the system is bistable, with an unstable fixed point at 

x*~ ^(1-^1 -4/3) + 0(a) (3.3) 

and two stable fixed points at 

x-^c7 + 0{a^), x+^ ^(l + v/l-4/3)+0(a) (3.4) 

This is the regime of interest as we wish to characterize the transition times 
between the two stable fixed points when the system is stochastic with weak 
fluctuations. 




Figure 1: Bifurcation diagram for the deterministic dynamics. 

3.2 Quasi-deterministic limit Og — )■ oo 

A velocity jump process, as described in the Introduction, can be obtained by 
taking the limit ac — > oo in either the discrete or semi-continuous process (both 
processes converge in this limit). This limit is discussed in [17] and later a 
metastable analysis was introduced in [28]. In this limit, the differential Kol- 
mogorov equation for both the discrete and semi-continuous processes converges 
to 



d_ 
di 
d_ 
di 



Po(a;,i) 
pi(a;,i) 



d_ 
dx 
d_ 
dx 



[{a - a;)po] - a-, (x^po - /3pi) 
[(1 - a;)pi] + a; (x^po - /3pi) . 



(3.5a) 
(3.5b) 



It is worth pointing out that the QSS approximation converges to a very different 
process in this limit; the differential Kolmogorov (in this case Fokker-Planck) 
equation becomes 



du 



-^((v(.)).) + -- 



D{x) 



du 
dx 



(3.6) 



where (v(a;)) and D{x) are given by (2.19) and (2.21), respectively. Moreover, 
one can derive (3.6) as a QSS approximation of (3.5) using the same projection 
methods as were used to derive (2.18). 



3.3 Adiabatic limit aj — )• oo 

The reduction procedure is based on a projection method very similar to the 
QSS reduction in Section 2.2. We leave the details to Appendix B and state the 
result. The limiting master equation is 



du 

'dt 



pKu, 



(3.7) 



10 



where 

1 
u{n,t)=^Ps(n,t), (3.8 

s=0 

and 






^(E+-l)f + (E--l)/(-;^), f{x)^"-i^^. (3.9) 



Note that at deterministic fixed points, x^ c = *,±, /(a^c) = a;^. 

4 Long-time asymptotic approximation of the prob- 
ability density function 

So far we have introduced the model problem, the discrete process, along with 
two approximations, the semi-continuous and QSS processes, which account 
for both noise sources. We also presented limiting processes that result when 
one of the two noise sources is removed, resulting in either the adiabatic or 
QD limiting processes. We now present a systematic perturbation method to 
analyze metastable, or long-time, behavior of the discrete and semi-continuous 
processes. 

Suppose we have a CK equation of the form 

d 

— p(a;,s,i) = -£eP, (4.1) 

where C^ is a linear operator acting on state variables (x, s). The state variables 
are split so that cc € r2 C K represents the external state and s € N represents 
the internal state. Note that one can easily generalize this theory to the case 
a: G fi C 'MJ^ (see [28]). For illustration, take C^ to have the form 

d d"^ 

£eP = -Q-{v{x, s)p) - £^(^(2;, s)p) 

^ ("(^' ^l«')p(a;, s', t) - a{x, s'|s)p(x, s, t)) , (4.2) 

s' 

where e ^ 1 is a small parameter (we will define e for the model problem later 
in this section). The CK can also be written in conservation form with 

-u{x,t) ^ - — J{x,t), (4.3) 

where 

u{x, t) = Y^ pix, s, i), J{x, t)^^i v{x, s)j>{x, s, t) - e— 6(a;, s)p(x, s,t)\ . 



11 



(4.4) 

Assume that C^ has a complete set of eigenfunctions, {0j(a;, s)} and adjoint 
eigenfunctions {^j(a;,s)}. If the initial condition is p(a;, s, 0) ~ (5(a; — a:;o)(5(s — sq), 
the solution can be written 



p(x,s,i) = V^j(a;o,so)(/)j(x,s)e ^^■*, (4.5) 



where we assume that all of the eigenvalues, Aj, are nonnegative. Since we are 
interested in metastable behavior, assume that in the limit e ^ 0, a determin- 
istic, bistable process is reached for the external state x. Label the two stable 
fixed points x± and the unstable fixed point a;* and assume x_ < x* < a;+. 

The random process will look very different if the external state starts at 
xq < Xs, or cco > x.^,. For the sake of illustration assume that xq — x^. On inter- 
mediate time scales, the solution will converge to a stationary density around 
X- that, figuratively speaking, does not see the other stable fixed point — or said 
another way, the solution does not see beyond x* . Slowly, over a long timescale, 
the solution converges to the full stationary density as probability slowly leaks 
out past x* toward x^. The timescale for this long-time convergence is expo- 
nentially large (i.e., 0(e^/^)). Since a stationary solution exists, the smallest, or 
principle, eigenvalue is Aq = 0, and the stationary density is the eigenfunction 
0o(a;, s); that is, we normalize the principal eigenfunction so that ((^Oil) = 1 
with respect to the inner product defined by 

{I{x,s),g{x,s))= I ^f{x,s)g{x,s)dx. (4.6) 

Jxen g 

The separation of time scales in the problem can be exploited to approxi- 
mate the solution. To understand how this works consider the process where a 
boundary condition is placed at x* so that the process truly does not see beyond 
the unstable fixed point. We want to consider two different boundary condi- 
tions: refiecting and absorbing. To distinguish between each case, we write the 
principal eigenvalue and eigenfunction as A*^"-*, 0^"-' and A'^''-', (j)^^"^ for absorbing 
and refiecting boundary conditions, respectively. If we place a refiecting bound- 
ary at X* the principal eigenvalue A^'"' = 0, but the eigenfunction (j)^"^' is now 
restricted to x G f7_ = (— oo,x,t) (or x G 51+ = (x,t,oo) if we instead assume 
that Xq > x.^,). We call this the quasi-stationary density. 

Now suppose that an absorbing boundary is imposed at x*. In this case, 
no stationary density exists, and the principal eigenvalue is perturbed by an 
exponentially small amount, that is, A*^") = 0(e^'-^/'^), for some C > 0. The 
eigenfunction (/)(''^ is also perturbed, but away from the boundary, (^^"^ ~ (j^^^K 
which turns out to be straight forward to compute using a WKB approximation 
method. Thus, if we can calculate the eigenvalue and eigenfunction, we have an 
accurate approximation to the absorbing boundary problem with 

p(x, s, t) ^ 0(") (x, 5)e-^'°'*, t\["-^ » 1, (4.7) 
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or, since <p^'^^ ^ 4>^'^\ 

p(x,s,i)-(^M(a;,s)e'^'°'*, iA^°^»l, e « 1. (4.8) 

We discuss how to approximate A*^"-* later in this section. 

This approximation can be repeated for the initial condition xq > x+, and 
a different principle eigenvalue and quasi-stationary density are obtained, call 
the eigenvalues Aj? and quasi-stationary densities ^^T (x, s). The full system, 
without any boundary condition imposed at x*, can then be approximated by 



1 ] q-{i)4>_{x,s), X < x^ 
q+{t)(j)y{x,s), X > x^ 



p(x,s,i)~^<i ,, (,) , , , (4.9) 



where q± (t) satisfy the system of ordinary differential equations 

^ = -Xi'\-+X^:\+ (4.10) 



dt 



= AL''^<Z_-Afg+, (4.11) 



with (?_(0) = 1 and q+(0) = if a;o < a;*, or g_(0) = and g+(0) = 1 if Xg > x*. 
To obtain an approximation of A^"' , we use a spectral projection method that 
makes use of the adjoint operator £* . This method was first developed for scalar- 
value PDE eigenvalue problems [15,21] and later generalized to a vector-valued 
PDE eigenvalue problem [16,28,30]. The present treatment further generalizes 
the method. Consider the adjoint eigenfunctions {^j}, j = 0, 1, • • • , satisfying 

/::0 = A,^„ (4.12) 

and take {(f>i,^j) — Sij so that the two sets of eigenfunctions are biorthogonal. 
We use the same notation to distinguish between the two boundary conditions 
for the adjoint eigenfunction. If the boundary is reflecting, the first adjoint 
eigenfunction is <^'^'") = 1, and if the boundary is absorbing then ^*^°) ^ ^'^''^ away 
from the boundary, but develops a boundary layer at x* . Using integration by 
parts we have 

(J (/)('■), A^")^'"^)^ = (J^W,/::^^")^ = ^£,0M,e*'^^) + J(0^''\e^"^), (4.13) 
where the boundary contribution, 

J(0W,e('^)) = e^6(2^„s)</)M(x„s)^e^'^H^*,s)- (4.14) 

is nonzero because 0''"' does not satisfy the absorbing boundary condition. 
Then, since £^4''^'^^ = 0, the principal eigenvalue is 
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The above identity can be used to approximate the principal eigenvalue as fol- 
lows. Since away from the boundary x*, ^^°-' ~ ^^''^ = 1, we can make this 
substitution for the term in the denominator of (4.15) so that 



Xi-) 



Ji 



(4.16) 



with exponentially small error. Notice that the denominator is then well ap- 
proximated by the normalization factor for the eigenfunction <j)'^^\ We cannot 
make the same substitution for the term in the numerator, since (4.15) becomes 
a formula for A*^*"^ — instead of A^°^. Of course, in some sense zero is actually 
a very good approximation because the error is 0{e~'^^'^), but to capture the 
metastable behavior we need to capture the small exponential. Note that we 
could have just as well used (j)^""^ and ^'^''■' in (4.13) instead of ^^''^ and ^('''. We 
choose the later because it simplifies the boundary layer analysis. 

The recipe for approximating the solution requires approximations of the first 
eigenfunction and the first adjoint eigenfunction, where the latter satisfies the 
appropriate adjoint absorbing boundary condition. For the examples considered 
here, the eigenfunctions are treated as vector-valued functions of x. That is, the 
internal state is restricted to s £ {0, 1} and we write (j)o{x,s) = [(/>o(.t)]s with 
^0 :K^^M2. 



4.1 WKB approximation of the eigenfunction 



Or), 



X] 



Before proceeding with the perturbation analysis, it is convenient to make a 
few definitions. The CK equation for either the discrete process (2.4) or the 
continuous process (2.16) can be written in matrix form. Specifically, we are 
interested in the properties of the transition rate matrix for the internal state. 
Based on the transition rate diagram (2.13) the transition rate matrix is 



Aix) 






(4.17) 



Transition rate matrices like this are members of a family of matrices called 
W-matrices. Notice that the columns sum to zero, which means that the matrix 
is singular and the vector, 



(4.18) 



is the left eigenvector corresponding to a zero eigenvalue. For a transition rate 
matrix to be a W-matrix, it must have negative diagonal elements, nonnegative 
off-diagonal elements, and it must be irreducible. One can show, using the 
Perron-Frohenius theorem, that the nuUspace of a W-matrix is one dimensional 
and that the right nuUvector has strictly positive elements. We denote this 
vector as p, where Ap — 0. Note that if the matrix ^ is a function of x then so 
is p. For the matrix (4.17) we have that 



p{x) 



/3 



(4.19) 
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where we have normalized p so that its elements sum to unity. For a simple 
two-state Markov process, where the external state is fixed, p is the steady state 
distribution. 

The WKB analysis is carried out in the small variable e defined as follows. 
Assume that the ratio between the two large parameters, a^ and aj, is fixed so 
that a-i = Lpac for fixed ip and define the small parameter e — l/ai — l/(ipac). 
It is convenient to put the CK equation (2.4) or (2.16) in matrix form so that 
the eigenfunction (up to terms exponentially small in e) satisfies 



[A(.T) + Sd]0o(a;)-O, 



(4.20) 



where A, given by (4.17), is the transition rate matrix for the internal state and 
d is a vector of linear operators governing transitions between external states. 
For the discrete process we have 



ddi 



1 



disc 



(e^^- l)x + ais-^"" ~ 1) 
(®^=" - l)x + {<B-^'' - 1) 



where the jump operators, 

e±^^/(x)^f:^/(")(.) = /(.±^.). 



(4.21) 



(4.22) 



r!=0 



are defined in terms of a Taylor series expansion. For the semi-continuous 
process 



d.s, 









where the drift and diffusivities are contained in 



v(x) 



a — X 
l-x 



, Hx) 



{(T + X)ip 

2 

{l+x)ip 

2 



We assume that the eigenfunction has the following WKB form 

1 



(t>o{x) ^ (ro(a;) + eri(a::) H ) exp 



-<i>{x) 



(4.23) 



(4.24) 



(4.25) 



where $ is a scalar functions and rp^i are 2- vectors (with Tq positive). Substi- 
tuting (4.25) into (4.20) and collecting leading order terms yields 



0(1) : [A{x) + Sh(.,*'(.))] ro{x) = 0, 
where, for the discrete process 



(4.26) 



hdisc(a;,p) = -(e^P - 1) (v(0) - xe'^P) , 



(4.27) 
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and for the semi-continuous process 
Kc{x,p) =pv{x)+p'^h{x). 



(4.28) 



with p = $'. An equation for $' is obtained by taking the determinant of the 
matrix in (4.26). We define the Hamiltonian^ as 



H{x,p) = hi(x,p)h2{x,p) - I3hi{x,p) ~ x^h2{x,p). 



(4.29) 



The equation for $' can then be expressed as H{x,^'{x)) = 0. In particular, 
for the discrete process we have 



Hdisc(x,p) = ^(e-'^P - 1)' + ^(e^P - ly 

--(/3a + x2)(e^P-l), 

and for the semi-continuous process 

Hscix,p) = bi{x)b2ix)p'^ + {bi{x)v2{x) + b2{x)vi{x))p^ 

+ {vi{x)v2{x) - i^biix) + x^{x)b2{x)))p'^ 

-{(3vi{x) + x'^V2ix))p. 



(4.30) 



(4.31) 



All of the roots of Hsc — are real, but there is only one nontrivial root that 
vanishes at the deterministic fixed points, namely 



$'(x) = 



-1 



3bi{x)b2{x) 



b,{x)v2{x)+b2{x)v,{x) + U\Zix)\'/'+-^^ 



X f cos(-arg(Z(a;))) - \/3sin(- arg(Z(a;))) 



where 



Z{x) = \ (C2{X) + 1^ \Cr{xY ~ C2{xf) , 



, (4.32) 



(4.33) 



Cx{x) = {bx{x)v2{x) + b2{x)v^{x)f 

-Zb^{x)b2{x)\v^{x)v2{x)~{x%2{x)+iib^{x))\ (4.34) 



^So named because of the similarities to classical Hamiltonian dynamics, even though it 
does not necessarily represent a physical energy. It does, however, define a scalar potential <&. 
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and 



C2{x) = -2{bi{x)v2{x) + b2{x)vi{x)f 

+ 9bi{x)b2{x){bi{x)v2ix) + b2ix)vi{x)) 

X [ui(a;)w2(a;) - {x'^b2ix) + l3bi{x))] 
+ 27{bi{x)b2{x))\x\2{x) + Pvi{x)). 



(4.35) 



For the discrete problem, the Hamiltonian "Hdiscla^jP) can be transformed to 
a cubic polynomial in g = e^'^. Then, the solutions are given by the positive 
real roots of the cubic polynomial 



(9-1)' 



f 



+ x{l + x + Lp{(3 + x^)+a)q-x^. (4.36) 
There is a single suitable root, given by 

Qi{x) \ 



q{x) 



-1 
3^ 



a;(l + Lpx) — a{l + x + Lpj3) 



1 



\Y{x)\ 



1/3 



\Y{x)\ 



1/3 



X ('cos(i arg(y(a;))) - V3sin(i arg(r(x)))') 



, (4.37) 



where 



^(a;) = \ {Q2{x) + i^^Qi{xf-Q2{xY) , (4.38) 

Qi{x) = (a::(l + -^x) + cr(l + .t + tpl3)f - 3crx(l + a; + (^(/3 + x^) + a), 



(4.39) 



and 



p(a::) = 9cra; (ct (1 + a: + /3(/3) + a; (1 + (y9x)) (l + cr + a; + (y3(;9 + x^)) 



- 2{cT [l + X + Pif) + X {I + Lpx)Y - 27cr^a;^ (4.40) 

The potential function <^(x) for both approximations is obtained numerically 
by quadrature. We rewrite the remaining term in (4.25) as 



ro(a;) = fc(a;)w(x), 



(4.41) 



where the approximation for the conditional internal state distribution (2.7) is 
determined by calculating the nuUspace of the matrix [A + Sh), which is 



w(a;) 



-h2{x,^'{x)) 

hx(x ,^' (x)) — h2(x ,^' (x)) 

h-i{x,^' (x)) 
hi(x ,^' (x)) — h2{x,^' (x)) , 



(4.42) 



The scalar function k{x) is a normalization factor and is determined at higher 
order. 
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To calculate k{x), substitute (4.25) into (4.20) and collect 0(e) terms to get 



0(e) : [A + Sh(2;,$'(x))] ri = -pSh^w+fc 

The left nuUvector, satisfying i-'"[A + T^h(x,<s>'(x))] = 0, is 

' h2(x,^'{x)) ' 

l{x) = 1 



Sh +-^"(x)T.h Iw + Sh^ 
tip:. I 2 V y Upp I Up ^_^ 



h2{x,^'(x)) 

l3+x^ 
hi{x,^'(x)) 

P+x^ 



(4.43) 



and it follows from the Fredholm Alternative theorem that k{x) satisfies 
dk 



dx 



+ ^'{x)k = 0, 



(4.44) 



where 



^'(x) = 



F{x)Up,ix, ^'(x)) + '^<S?"ix)F{x)Upp{x, ^'{x)) 



(4.45) 



F(a;)Hp(a;,$'(a;)) 
for X ^ x±,x^, and 

H(x,p) = [A + Eh(.,p)]w(x). (4.46) 

We can express $" (x) in terms of partial derivatives of the Hamiltonian with 



Hx{x,^'{x)) 



(4.47) 



For more about evaluating the limit x — >■ Xc, Xc — x±,x^,, of $"(a;) and ^'(x), 
see Appendix C. Integrating (4.44) gives 



k{x) — exp [— $(a;)] , 



(4.48) 



where \E'(a;) is computed by numerical quadrature of (4.45). Finally, to leading 
order, the eigenfunction is 



4>o{x) ^ w(a;)fc(x) exp 



--$(x) 



(4.49) 



where k{x) and \v(x) are defined by (4.48) and (4.42), respectively, while $(a;) 
is determined by numerical quadrature^ of $'(x), defined by (4.32) for the semi- 
continuous process and (4.37) for the discrete process. Note that at fixed points, 
Xc = x±,x^,, we have that w(a;c) — p{xc), where p{x) is defined by (4.19). 

^In practice, we find that the best way of numerically integrating <I>'(x) and '^'{x) is to use 
Chebychev approximation methods (we use the GNU Scientific Library). 
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4.1.1 Adiabatic limit 

In the adiabatic limit, we have that w{x) = p{x) (see (4.19)). The quasi- 
stationary analysis is well known for the reduced process, and because of the 
reduction in the complexity of the process, a simple analytical approximation 
can be obtained. We quote the result here and refer the reader to [5,32]; that 
is, 



and 

k{x) = ^ , (4.51) 

Va;/(x) 

where f{x) is given by (3.9). 

4.1.2 Quasi-deterministic limit 

Both the discrete and semi-continuous process converge to the QD process (Sec- 
tion 3.2) in the limit 99 — > (i.e., ac ^- 00 with a; fixed). A similar analysis 
can be carried out on the CK equation (3.5) for this process (see [16,28,30] for 
details). The result is a simple analytical approximation. For cr < a; < 1 we 
have 



ri- 


-x~\ 


1- 

X- 


-a 
-a 


Li- 


-(7- 



w{x)^ Iz- , (4.52) 



$(a;) = -/31n(l ~ x) ~ a^ \n{x - a) - -{x - af - 2a{x - a), (4.53) 



and 



k{x) = \- -. (4.54) 

[x ~ (J){1 — X) 

4.2 Singular perturbation approximation of the adjoint 
eigenfunction ^*^"^(x) 

The primary purpose of this paper is to develop the singular perturbation tech- 
niques to approximate the adjoint eigenfunction. The WKB method used in 
the previous section is a straightforward extension of those presented in ]28] for 
the QD limit, and it provides only an approximation of the stationary density, 
not the timescale for metastable transitions (i.e., the principal eigenvalue A). 
To get information about transition times we must have an approximation of 
the adjoint eigenfunction. We emphasize that singular perturbation methods 
are ideally suited to higher dimensional problems that display nonequilibrium 
stationary behavior. 
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To make the following analysis applicable to a more general class of Markov 
processes, we assume that there are M internal states; for the discrete and semi- 
continuous process M — 2. Here we assume that h G M.^ , and the transition- 
rate matrix ^ is M x M. Because the boundary layer analysis of the semi- 
continuous process is very different from the discrete process, we treat the two 
cases separately. 

4.2.1 Semi-continuous process 

Up to terms exponentially small in e, the first adjoint eigenfunction satisfies 



[A{x)]^ + eS^(j,)— + e^Sb(a;) 



'^ '" '^' U^) = 0, (4.55) 



dx dx'^ 

along with the absorbing boundary condition, 

^o(O) - 0. (4.56) 

The outer solution, which does not satisfy the boundary condition, is exactly 

Cout = 1- 

To obtain an approximate solution that also satisfies boundary conditions, 
we must rescale x ^ x^ -\- e^ z, for some 9 > 0. A reasonable first try is to take 

= 1 so that X — x^ + ez. Equation (4.55) becomes 

[^(x*)]^ + Ev(..)^+Sb(..)^ ^bi(z)=0. (4.57) 

The solution is a linear combination of the subsolutions CjT je~^^^ , j = 0, • • • , 2M- 

1 where 

{A{x,f - 7,Ev(..) + lpb{x.)) Tj = 0, (4.58) 

and Cj, i — 0, • • • , 2Af — 1 are unknown constants. We can specify the first 
solution as To = 1 and 70 = 0. We also require the solution to be bounded 
in the limit z — )■ cxd, which means that Cj = if 7j < 0; although we do not 
know apriori how many of the eigenvalues are negative. Note that the boundary 
condition (4.56) provides a system of M linear equations for the 2M unknowns, 
Cj, which means that constraints to eliminate the remaining M unknowns are 
required to close the system. One such constraint eliminates an unknown (i.e., 
Co) by matching to the outer solution, leaving M — 1 more constraints we must 
find, which suggests that there are M — 1 negative eigenvalues. For simplicity, 
we order the eigenvalues so that 7^ < for j = M -|- 1, • • • , 2M — 1. 

For the moment, consider the matrices in (4.58) as depending on x so that 
Tj and 7^ are also functions of x. It is simple to show that To(a;) — 1 and 
7o(a;) = even if the matrices are evaluated away from x*. However, one of the 
solutions, label it j = 1, is 71 = p{x) — ^'{x) from (4.32). Moreover, Ti(x) — >■ 1 
and 7i(x) — > as a; — > a;*. In fact, 71 (x) vanishes at all of the deterministic 
fixed points since '^' [xc) = 0, for x = x±,x.^,, and is nonzero at all other points. 
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It follows that the zero eigenvalue has a degenerate eigenspace, and the solution 
must include a secular term involving the generalized eigenvector satisfying 



Mx*)'^C = Sv(a;.)l = v(a;*) 



(4.59) 



One can show [30] that the deterministic fixed points are the only points where 
this matrix has a degenerate eigenspace associated with the zero eigenvalue. 
The solution to (4.57) is thus 



M 



^bi(z) = cqI + ci{C - zl) + ^CjTjC 

i=2 



-7j2 



(4.60) 



However, because of the secular term, the solution remains unbounded in the 
limit z — >■ oo, and as a result, it cannot be matched to the outer solution. 
Therefore, there is a transition layer which sits between the boundary layer and 
the outer region. 

To find the scaling for this transition layer, we change variables to a; = e^s, 

for < 6* < 1, and define ^e{s) = ^o(^)- Introduce the asymptotic expansion 



Equation (4.55) becomes 



{A{x,) + e''sA'{x,) + ---y 



(4.61) 



+ ei-''(Ev(,.)+e^sSv'(,.) + ---) 



ds 



d^ 
ds'^ 



0(1): A(x,)^|f(.) = 0, 



X (C^ (s) + e-i'-p (s) + e^'^if (s)) = 0. (4.62) 
Setting e = in (4.62) yields 

(4.63) 
which implies that 

if\s) ^ a^{s)l, (4.64) 

for some scalar function ao(s). The leading-order terms in (4.62) then become 

+ 0(e) + 0(e^+") + 0(e2(i-s)) + Oie^") = 0. (4.65) 
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Note that £^A^1 = for all n > 0. Setting k = 1 - 6* we have 

0{e^-') : A(x.f^i'\s) = -a^(s)v(x.), (4.66) 

and smce p{x^,)'^\'{x^,) — 0, the solution is 

il'\s) = -a',isK, (4.67) 

where ^ satisfies (4.59). Thus, 

$g{s)^ao(s)l-e'-'a'a{s)C. (4.68) 

The function ao(s) is determined at higher order; we find 

ei-«A(a;,)^|f ' + e'a'^ {z^r'{x,) - zA'{x,fC) 

- e^-'al (Sv(x.)C - b(x,)) = 0. (4.69) 

Setting e = 1/2 yields 

0{e) : A{x,fif\s) = a;;(s)S,(,.)C-aJ,(s)s (v'(x.) - A'{x,fC,) , (4.70) 

and the resulting solvability condition is 

"( \ ( {p{x^Yv{x^))' \ ,. . „ .. „,. 

Qo(g)-g I .T (s^ 7 — , / xN «o(g) =0- (4-71) 

\p{x.^y (Sv(^.)C-b(a:;,))y 

Note that Ap = Q^ A'p = -Ap'. Furthermore, C^A(x*)p'(a:*) = v(a;*)^p'(a:*). 
Hence, p(a;H,)'^v'(a;*) — p{x^)'^A'{x^)C = {p{x^)'^v{x^)y. One can show that 
(see Appendix A) 

^("*)'^"'("*) - -$"(..). (4.72) 



p{x.^.)'^{h{x.^.) - ^v(x.)C) 

Assuming that $"(a;*) < 0, the solution to (4.71) is 

a[,(s)=cie5*"(-*)«', (4.73) 

and 

ao(s)=co + ci r e"-'^"^^'^y"dy, (4.74) 

Jo 

where co,i are unknowns constants. Note that since x = a;* is a local maxima 
of $(a:), we assume that $"(a;*) < so that ao(s) — >■ as s — )• cx). The solution 
(4.68) becomes 

li/2(s) ^ (co + ci f^ e^*"(-)^^dy) 1 - 6V2cie^*"(-)^^C, (4.75) 
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which replaces the first two terms in (4.60) (i.e., cqI + Ci(C — -^1))- Notice that 
the solution is now bounded in the limit z — >■ oo, which allows us to match it to 
the outer solution; we require lims_>.oo ^1/2(2) = 1 so that 



Co + ci / e 
/o 

and take 

Co = 1 - ci 

As s -> 0, 



5*"(^*)'''ds 



Co + ci 



2|$"(x, 



1, 



2\^"ix,)\' 



li/2(s)-- (co + cis)l-e^/^ciC 



(4.76) 



(4.77) 



(4.78) 



This matches with (4.60) if Cq = Co and ci = — e^^^Ci. A uniform solution, valid 
throughout the boundary layer and transition regions, is 



U^) 



1-ci 



2|$"(x*)| 



{x-x,)/e^ 



,5*"(:r.)y= 



dy 



M 



J=2 



~lj{3:-x,)/e 



(4.79) 



The remaining unknown constants Cj, j — 1, ■ • • , M, are determined using the 
absorbing boundary condition (4.56), resulting in the linear system of equations. 



"^'ijwk:fr'"'')-p''''^'^ 



(4.80) 



Using this general result, we can calculate the adjoint eigenfunction for the 
model problem. We have that M = 2 for the semi-continuous process (Section 
2.1), and the solution to (4.80) is 



,.,/^^::M+o(eV2), 



TT 



C2 



-«"'>/^^(Aj-°'* 



where A = [1, —1]^. A simple calculation shows that 



c = 



1-x^ 



(3 + 1 



(4.81) 
(4.82) 

(4.83) 
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and 



To = 



P + xl 



-/? + 72(l-x0 + f7l(l + a^*y 



where T2 is normalized so that 

l2^l^b(x,)P{x*) = T^Ev(2:.)P(a;*) = 1. 

4.2.2 Discrete process 

The adjoint eigenfunction for the discrete process satisfies 

([^(x)]^ + I]d.)6(-^) = 0, 
where A is given by (4.17); the adjoint operator is 



(4.84) 
(4.85) 

(4.86) 



x(®-^"^ - 1) + o-(®^^"^ - 1) 
x(e-^^ - 1) + ((E^'' - 1) 



(4.87) 



(4.88) 



where <b ''^ are defined by (4.22). The absorbing boundary condition is 

6(n*) = 0, (4.89) 

Once again, the outer solution is ^out = 1- 

Motivated by the boundary layer analysis of the previous section, we rescale 
with a; = a;* + e^s. We are interested in two cases: 9 = 1 and 6 = 1/2. In the 
former case, the scaling simply returns the process to a discrete variable since 
X = ipen. Let fi = n ~ n^ and ^{,(n) = ^0(2;* + (pen). Then to leading order 

^[A{x,)f^hin)+x, (6(ri - 1) - Un)) + ^.io) (4(n + 1) - 6(n)) = 0. (4.90) 

Solutions have the form 

6(n) = T.fif. (4.91) 

Substituting this into (4.90) yields 

[ipii,[Aix,)f + (,,{^1, - l)E,(o) + x4fi, - 1)1] T, = 0, (4.92) 

with fij determined by the characteristic equation. 



det [(/7/ij[A(a;*)]^ + a*j(Mj ^ l)Sv(o) + a^*(Mi - 1)^] = 0. 



(4.93) 



Just as in the previous section, one of the linearly independent solutions is 

Un) = C- ^nl, (4.94) 

where ^ is given by (4.59). 
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On the other hand, if 6* = 1/2 we recover (4.62), which means that we can 
use (4.75) as a solution instead of (4.94). The uniformly valid solution is then 



U^) 



1-ci 



2|<i>"(x*)| Jo 



,^<S>"ix,)v' 



dy 



l/2^^gl*"(:..)(2i^fi^ 



M 



C + Y.c,T,^Ji^-' , (4.95) 



3=2 



where we assume that \fij\ < 1, j = 2, • • • , M. The unknown constants are 
determined by the boundary condition ^0(2;*) = 0, which results in the same 
linear system (4.80) as obtained for the semi-continuous process, except that T^- 
now satisfies (4.92). 

Turning to the model problem, calculation of ^0(2;*) for the discrete process 
is very similar to the semi-continuous process. Using (4.81), (4.82) and (4.83) 
in (4.95), we need only calculate T2 and ^2 satisfying (4.93) and (4.93). We 
find that 



P + xl 




-Vxl^i2 




xlR 


.(^A^i - 


- {a + x^)fi2 + x^) - 


- (Pxtn2_ 



and 



T2 = 



R — {afil — (cr + x^)^2 + a;*) — (p{/3a + x1)^2 



M2 



-Lp{l - a)xl - (cr -I- l)(x* - (j)x^ + x^\/U 



(4.96) 
(4.97) 

(4.98) 



2cr(o- — a;*) 
U = {a- l)(^^(cr - l)xl - 2ip{a + l)xl) 

+ ia- 1)((2ctV + (1 + 2ip)a - l)xl - 2a{a - l)x, + a^{a - 1)). 

(4.99) 

4.3 Principal eigenvalue 

Now that we have approximations for the right and left eigenfunction, we can 
construct the approximation of the principal eigenvalue using the spectral pro- 
jection method (see (4.15)) outlined in the introduction of this section. For 
brevity, we only show steps for computing the eigenvalue corresponding to the 
right well for which the minimum is x+. The corresponding result for the left 
well proceeds similarly, and the final formula differs only by substitution of x_ 
for x_|_. 

First, for the semi-continuous process, substituting the eigenfunctions (4.49) 



25 



and (4.79) into (4.14) yields 



2e|$"(a;, 



X k(x^) exp 



h{x^)'^ p{x*) + i ^Tj ] 72T^Sb(:r,)P(a:0 



-$(x* 



(4.100) 



2e|$"(a;* 



[(b(x*)) + D{x^.)] fc(a:*) exp 



;$(a;*) 



where p, C, and T2 are given by (4.19), (4.83), and (4.84), respectively. Notice 
that 






'y2'^2^b{x,)P{x*) 



{x,-a){l-x,) _ 

P + xl -''^''*^' 



(4.101) 



the diffusivity from the QSS reduction (2.21) evaluated at the unstable fixed 
point. (Clearly, the boundary term J{(f>o,^o) is related to the classical diffusive 
flux J{u) — D-j^.) The normalization constant is well approximated by 



Finally, the eigenvalue approximation is 

A.^(|^l*"(..)l*"(x,))|Mexp 



■(«(!.)- «(l±)) 



where 



B = (fx^, + 



(x* — cr)(l — a;*) 



/3 + x2 
The preexponential factor is equivalently 



k{x±) 



exp[-(*(x,)-*(a;±))]. 



(4.102) 



(4.103) 



(4.104) 



(4.105) 



Recall that $(a;) and 'i'{x) are determined by numerical integration of (4.32) 
and (4.45), respectively. 

The discrete version of (4.15) can be obtained using a summation by parts 
argument. The resulting boundary contribution is 

J(0o, 6) = Mx*f^.{oMx* + IK). (4.106) 

The first two terms in ^0(2^* + l/'^c) (see (4.95)) can be expanded in l/a^ <^ 1 
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(the third term is the boundary layer solution), and the result is 



(4.107) 



e^/'c. 






ipx^ + (/i2 - 1) ( ^^^ ) T^Sv(o)P(a:^*) 



X k{x■^,) exp 



-$(x* 



where T2 is given by (4.96). Similar to (4.101), we find that 



(M2 - 1) i ^T-f ) '^2^A0)P{x*) = D 



(x*). 



(4.108) 



The eigenvalue for the discrete process has the form (4.103) with the following 
substitutions. The potential function ^{x) is determined by numerical integra- 
tion of $'(x) = -"y^feU.^ where q{x) is given by (4.37), and \E'(a;) is determined 
by integration of (4.45), with h{x,p) given by (4.27). 

5 Results 

In this section, we compare the approximations of the stability landscape, de- 
fined as — eln(M(a;)), and of the mean time of a metastable transition from the 
minimum of one well to another. The shape of the stability landscape can be 
described as a double-well potential, and in Fig. 2 it is shown for a — 0.015 and 
two different values of the bifurcation parameter, /3, located within the region of 
deterministic bistability (see Section 3.1). The stability landscape is shown in 
two columns of plots, each using different parameter values. In the left column 
/3 = 0.24, which is near the bifurcation point that eliminates the right stability 
well, and in the right column /3 ~ 0.11, which is near the bifurcation eliminating 
the left stability well. Each row shows a different value of e with ip = 1 so that 
both noise sources are present. Approximations of the stability landscape are 
given by ^{x) + €'^{x), where <& and ^ are defined in Section 4.1. Note that 
the WKB approximation of the discrete process brakes down as x — > due 
to small copy number, requiring a boundary correction (see Appendix E). Each 
approximation — the QSA discrete and semi-continuous approximations, and the 
diffusion approximation — is compared to a numerical approximation obtained 
by SVD decomposition in the top two rows for which < e <C 1. In the bottom 
row we take the limit e — !• 0. Note that the SVD approximation cannot be 
computed for this case. First, we observe that the QSA approximation of the 
discrete and semi-continuous process are so close that they are indistinguishable 
for every parameter set. (Indeed, we find this to be the case for all of the re- 
sults presented in this paper). On the other hand, the diffusion approximation 
shows significant inaccuracies, particularly in the left stability well. The most 
significant aspect of the stability landscape that affects metastable transitions 
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Figure 2: The stability landscape, — elog(u(a;)); for Lp = \ and a = 0.015. 
Each row shows results for a different value of e, and each column shows a 
different value of the bifurcation parameter /?. The light blue curve is the 
quasi-steady state approximation, the green curve is the semi-continuous QSA 
approximation, and the red curve (which cannot be seen beneath the green curve 
because both approximations are very close) is the discrete QSA approximation. 
For nonzero e (the top four panes), the first and second order approximation 
$(x) -I- 0^{x) is compared to the value of — elog(ps) obtained by a numerical 
SVD decomposition, shown as 'x' symbols. For e = (the bottom two panes), 
the leading order approximation $(x) is shown. Note that the SVD solution 
can not be computed in the e — >■ limit. 
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is the height of each well in the e — ^ limit. Although the diffusion approxi- 
mation does show some error in right stability well, including the height when 
e = 0, these differences are much less significant than the differences in the left 
well region. Even for the left well, the diffusion approximation is not always 
inaccurate. Indeed, all of the approximations closely agree when e = 0.005 and 
P = 0.24 (first column, second row of Fig. 2). However, for other values of e 
(top and bottom row) this is clearly not the case. 

To examine the differences in the approximations more closely, we plot the 
absolute error in the stability landscape and the error in the conditional internal 
state distribution in Fig. 3 for the parameter values used in the left column of 
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Figure 3: Absolute error. A comparison of each approximation to the numerical 
SVD result, for /3 = 0.11, a — 0.015 shown in Fig. 2 (right column). The top row 
shows the error in stability landscape —e\og{u{x)) for e = 0.01 and e = 0.005. 
The bottom row shows the error in the internal state distribution for the same 
values of e. The colors for each curve are the same as in Fig. 2. 

Fig. 2 (i.e., a — 0.015, tp — 1, and /3 = 0.11). The conditional internal state 
distribution w(x) is (4.42) for the discrete and semi-continuous QSA approxi- 
mations and (4.19) for the QSS approximation. These are again compared to a 
numerical approximation obtained using an SVD decomposition, and the error 
is measured using the 1-norm (i.e., X^Lo l[wsvd(a;)]s - [wapprox(a:)]s|)- The dis- 
crete and semi-continuous QSA approximations of the stability landscape show 
errors primarily in the left well region, while the QSS approximation also shows 
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some error in the right well. Interestingly, the conditional internal state distri- 
bution error is significant for the QSS approximation, peaking at 25% between 
x* and x+. We expect this error to be quite small near the deterministic fixed 
points, where all the approximations agree. We emphasize as one of the key 
results of this paper that away from fixed points, the conditional internal state 
distribution is not always close to the steady-state distribution as assumed in 
the QSS approximation method. This has been shown rigorously for velocity 
jump processes [30], for which the QD limit is an example. 

The approximation of the mean time for a metastable transition between 
wells is shown in Fig. 4. The mean escape time approximations, defined as 
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Figure 4: Mean exit time approximations compared to Monte-Carlo results, (a) 
Exit from the left well for the same parameters as used in the first column of 
Fig. 2. (b) Exit from the right well for the same parameters as used in the 
second column of Fig. 2. 

T ~ 1/A (see (4.103)), are compared to exact Monte-Carlo (MC) simulations 
for parameter values used in Fig. 2. The mean escape time is plotted as a 
function of 1/e because each approximation is a linear function of this quantity, 
with a slope determined by the height of the potential well in the e — >■ limit (see 
Fig. 2 bottom row). Escape from the left well (for (3 = 0.11, Fig. 2 left column) is 
shown on the left, where the discrete and semi-continuous QSA approximations 
are in good agreement with MC simulations. The three approximations converge 
near e = 0.005 consistent with Fig. 2 (first column, second row). 

A somewhat unexpected result is obtained for escape from the right well 
(corresponding to the right column of Fig. 2). All three approximations are very 
close, and the QSS approximation is actually more accurate for smaller values of 
1/e. The difference in the slope of each approximation is slight (see Fig. 2 right 
column, bottom row) and the error in the QSS approximation should grow as 
1/e — )• oo. We cannot offer a definitive explanation for the accuracy of the QSS 
approximation for escape from the right well. One explanation consistent with 
our results is that the QSS approximation is valid for larger values of x, which 
seems reasonable since it relies on fast transitions between internal states and the 
rate of transitioning from the inactive to the active internal state is proportional 
to x^. However, this is inconsistent with the error in the conditional internal 
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state distribution shown in Fig. 3 (bottom row), which is the key assumption 
underlying the QSS approximation. 

Finally, we compare the mean escape time in the adiabatic limit aj — )• oo 
and in the QD limit etc — > oo. In Fig. 5, the mean time for escape from the left 
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Figure 5: The mean exit time for escape from the left well to the right well, 
with j3 — 0.23, (J = 0.04. Three different approximations (solid curves) are 
compared to Monte-Carlo simulation results (symbols). The discrete (red) and 
semi-continuous (green) QSA approximations are indistinguishable. Also shown 
is the QSS approximation (light blue), (a) The mean exit time as a function 
of tto for fixed a-^ — 333. (b) The mean exit time as a function of a; for fixed 
ao = 200. 

well is shown for a = 0.04 and j3 = 0.23. In contrast to previous results, we 
do not fix (/5 = a-Jae — 1. Fig. 5 (right) illustrates that the discrete and semi- 
continuous approximation converge in the QD limit, and as expected, the QSS 
approximation error is significant. In the adiabatic limit (Fig. 5 left) the discrete 
and semi-continuous QSA approximations show close agreement for all values 
of a-i, and as expected, all three approximations converge as ai — > cx). Even 
though the discrete and semi-continuous QSA approximations do not converge 
in the adiabatic limit, the difference is very small. This suggests that a diffusion 
approximation for the external state — recall that we used such a procedure 
to derive the semi-continuous process from the full discrete process — may be 
valid in certain situations, which is interesting since diffusion approximations 
generally break down for metastable behavior due to large deviation errors. It is 
possible that the good agreement that we see for the example problem is due to 
the linear nature of the birth-death process governing transitions in the external 
state (i.e., that it is due to the simplicity of the example problem). Since the 
general QSA procedure presented here does not depend on this assumption, it 
would be interesting to see how this type of diffusion approximation behaves if 
the external state is governed by a more complicated process. 
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A Curvature prefactor 

The purpose of this section is to show that the part of the eigenvalue estimate 
that contains information about the curvature of the stabihty well at the stable 
and unstable fixed point is unaffected by the diffusion approximations in Sections 
2.1 and 2.2. This is a reflection of the fact that diffusion approximations, in 
general, are accurate in a neighborhood of deterministic fixed points. Each 
eigenvalue approximation (4.103) and (2.26) contains a prefactor term of the 
form •\/|$"(a;*)| $"(x±). We would like to show that, when evaluated at a fixed 
point, call it Xc, the second derivative of the potential function for the discrete, 
semi-continuous, and QSS processes are all identical. Clearly, evaluating the 
derivative in the first two cases is messy. However, we can also express the 
second derivative in terms of the Hamiltonian functions (4.53) and (4.30) as 
follows. 

Differentiating T-L{x, $'(a;)) = with respect to x yields 

^n{x, <^'{x)) = n^{x, $'(a:)) + <^"{x)np{x, <^'{x)) = 0, (A.l) 

and it follows that 

np(x,^'{x)) 

However, we have that 

■Hpixc, 0) = n^{xc, 0) = n^^ixc, 0) = 0. (a.3) 

A formula valid at fixed points can be obtained as follows. Differentiating 
H{x, <&'(a;)) = twice with respect to x yields 

^n{x, $'(a;)) = n^, + <P"n^p + $"(Hp, + '^'"Hpp) + <^""Hp - 0, (A.4) 
and it follows from (A.3) that 

cl,"(a;j = ^f- ^ " ^ (A.5) 

§,nix,,o) 

A straightforward calculation shows that Taylor expansions of both Hamiltonian 
functions (4.31) and (4.30) agree to 0{p^), and that 

which is also the second derivative of (2.25) since {v{xc)) = 0. 
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B Adiabatic limit of the discrete process 

Consider the Master equation for the probabihty distribution function 

Pj{n,t) = p{j,n,t\jo,no,to). 
In matrix/operator form, the CK equation is 

|=X..+ ^X.p, (B.l) 

where p(n,i) — (pi(n, t), p2(n, t), • • • , pM(n, t))-^; Li — Sp is a diagonal ma- 
trix of linear operators acting on n, each of which has a W-matrix representation; 
L2 is an Af X M W-matrix governing the transitions between internal states, 
with transition rates that may depend on n. Define the projection operator 

V = pl^, (B.2) 

where L2p{n) = 0, with p(n) > and J2i=i Pji'^) — li ^^i^ 1 = (1, 1, • • • , 1)^. 
We assume the solution has the following form 

p(n, t) = 7'p(n, t) + (/ - 7')p(n, t) = u{n, t)p{n) + ew(n, t), (B.3) 

where 

u(n,t) EE l^p(n,i), l^w(n,i) = 0. (B.4) 

Applying the projection operator to both sides of (B.l) yields 

-^p = VLi(up + ew). (B.5) 

at 

On the other hand, applying the orthogonal projection yields 

e-^ - e(/ - V)Liw ^ {I ~ V)Li{up) + L2W. (B.6) 

After setting e = in the above equation we get 

win,t) ^ -L^\I -r)L,{u{n,t)p{n)). (B.7) 

Substituting (B.7) into (B.5) yields the scalar-valued operator equation for 
u{n,t) 

^ = l^Liiup) ~ el^LiL^\l - V)Li{up). (B.8) 

One can rewrite (B.8) in matrix form to obtain a linear system of ODEs for the 
vector u(i) with elements Un{t) = u{n,t) 

ft - «"■ <■' ») 
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where 



M 



w^^E^. 



jpj- 



(B.IO) 



j=i 



Note that the operators 2?j are now exphcitly assumed to be matrices. The 
resulting transition rate matrix, W , reflects the reduced graph structure. The 
transition rates contained in W can be interpreted as the average of the tran- 
sition rates contained in Dj with respect to the steady state distribution pj. 
Note that since p > and T>j are W-matrices, it follows that ly is a W-matrix. 
In general, the reduced equation represents a Markov process only at leading 
order. 



C WKB/KM expansion 

Consider the action of the operator <b on g{x)e^'^'''^^^' where g{x) is scalar 
function and ^{x) ~ (f^{x). We have that 



^±dx 



g{x)e 



-ao#(x) 



E 

CO 



[gi^ 



(±1)" d" 



^-Qe*(rr 



(n-fe) 



(x) 



dx'^ 



,-ao*(a:) 



EStE : .9'"-'^Hx)(l + B.(-..$(..)))e--*(^), 



(±l)".5r- in' 

n=0 <= fe=0 

where Bfc is the A;th complete Bell polynomial. 



B„(/(a;)) = dot 



/' ("7')/" ("2')/*'^ ••• /'"^ ■ 

-1 /' ••• /("-2) 



/' 



(C.l) 



(C.2) 



and Br, — 0. One can show that 



^k{-aM^))^al{-^'{x)f-a',-'-{k-l)<^"{x){-'^'{x))'-^ + 0{a^-^). 

(C.3) 
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Expanding (C.l) in terms of 1/ao yields 



[9{x)e 



e±^^(q(x)e-"<^*(^) 



^ gT4'(^) 



g{x) ~ ^ U{x) T \9{x)'b'\x) ) + 0(a-2) 



c*(2;) 



(C.4) 



D Evaluating ^'{xc) for Xc = a:±,a:+ 

Using L'Hopital's rule, we find that 



*'(Xe) 



Hp^^ + -^"{xc){3Hpp^ + <i>"{xc)Hppp) + -<i>"'{xc)Hpp 



+l\xcfUp,{x,,0) + -^"ix,)l'{xcflipp{x,,0) 



i{xcfUp{xc, 0) + Hp,, + $"(xc)i/i 



pp 



(D.l) 

where li{x,p) is defined by (4.46) and partial derivatives oi H{x,p) = l'^ll{x,p) 
are evaluated ai x = Xc and p = 0, as for example, 



Hxp — 1 
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dxdp 



H(.Te,0). 



We also have that ^"(xc) is given by (A. 5), and 



$"'(Xc) 



^ npxxjxcO) + l<i>"{xc)nppp{xc,o) 

ttpp[Xc^ yJ) 



(D.2) 



(D.3) 



Note that H{x,p) ^ 'H{x,p), where 'Hlx.p) is the Hamiltonian (4.29). 



E a; — > limit of the quasi-stationary density 

The WKB approximation (4.49) of the discrete process brakes down in the 
limit x -^ Q, due to small copy number effects (i.e., fluctuations are on the same 
order) . This fact is not relevant if one is interested only in approximating the 
mean exit time. However, we also approximate the effective potential. Although 
4>(a;) is bounded in the limit a; — >■ 0, ^(x) has a logarithmic singularity. To 
correct this, we use the discrete master equation (2.4) to calculate ^o(O)- with 



</.o(0) = -(aiA(O) - aeI]v(o))"'</'o( — )• 



(E.l) 
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The WKB approximation can be use for the value of the eigenfunction near the 
boundary, with 



4>o{l/ac) = vf{l/ac)k{l/ac) cxp 



1 



$(iK 



Hence, 



(E.2) 



0o(O) = -{aiA{0) - q:oSv(o)) "^w(l/Q;e)A:(l/Q;o)exp 



--$(iK 

e 



(E.3) 
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